# importing the necessary libraries & Packages:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
import pypsa
import cartopy.feature as cfeature
import cartopy.crs as ccrs
import cartopy
import warnings
import atlite
import rasterio
from rasterio.plot import show
import xarray as xr
import netCDF4
warnings.filterwarnings('ignore')
from atlite.gis import ExclusionContainer
from atlite.gis import shape_availability
import rasterio as rio
import os
import country_converter as coco

import logging
logging.basicConfig(level=logging.INFO)

from shapely.geometry import Polygon, MultiPolygon, Point

# extracting GADM data belonging to Slovakia & creating the shape of the country including its regions:

svk_regions = gpd.read_file('gadm_410-levels-ADM_1-SVK.gpkg')
svk_regions = svk_regions.set_index('NAME_1')
display(svk_regions)
shape_svk = svk_regions.to_crs(3035).geometry     # Slovokia is in Europe (crs = 3035)
shape_svk.plot()

<Axes: >

###### !!! this section is performed only for presentation purposes !!   #########

df = pd.read_csv('global_power_plant_database.csv')     # Reading the data frame included all power plant data of all countries
df_svk = df[df.iloc[:, 0] == 'SVK']                     # extracting the Slovakia-related data
display(df_svk.head(5))
df_total_generation_groupby = df_svk.groupby('primary_fuel').estimated_generation_gwh_2017.sum()
#df_total_generation_groupby.head(5)

# Plotting the pie chart of total generation per technology in Slovakia:
plt.figure(figsize=(8, 5))
plt.pie(df_total_generation_groupby, labels=df_total_generation_groupby.index, autopct='%1.1f%%', startangle=90)
plt.legend(df_total_generation_groupby.index, bbox_to_anchor=(1, 0.5), loc="center left", title="Energy Carriers")
plt.title('Estimated Generation GWh (2017) by Primary Fuel for Slovakia')
plt.savefig('Estimated Generation(2017).png')
plt.show()
print()
print()

geometry = gpd.points_from_xy(df_svk["longitude"], df_svk["latitude"])
gdf_svk = gpd.GeoDataFrame(df_svk, geometry=geometry, crs=4326)

<GeoAxes: >

# Reading the file included all of the Exclusive economic zones of the countries, which is not the case for Slovakia because it is a landlocked country:
# But the files are read anyway:

all_EEZ = gpd.read_file('eez_v11.gpkg') 
all_EEZ = gpd.read_file('eez_boundaries_v11.gpkg')

## !!! For presentation purposes !!! ###

# Reading the Slovakia Shape from data "country_shapes.geojson" file:

url = 'country_shapes (4).geojson'
countries = gpd.read_file(url).set_index("name")
#countries.tail(5)
gebco_svk = "GEBCO_2014_2D-SK.nc"
gebco = rasterio.open(gebco_svk)
band = gebco.read(1)
fig, ax = plt.subplots(figsize=(8, 12))
countries.loc[["SK"]].plot(ax=ax, color="none")
show(band, transform=gebco.transform, cmap="RdBu_r", ax=ax)

<Axes: >

# Defining Exclusion and Inclusion areas as data frames:

Exclusion_area = pd.DataFrame({"land classes" :[111, 112, 121, 122, 123, 124, 142, 511, 512, 131, 322],
                   "code value" : [1, 2, 3, 4, 5, 6, 11, 40, 41, 7, 27],}, 
                    index = ['Continuous urban fabric', 'Discontinuous urban fabric', 
                  'Industrial or commercial units', 'Road and rail networks and associated land', 
                  'Port areas', 'Airports', 'Sport and leisure facilities', 'Water courses', 'Water bodies','Mineral extraction sites','Moors and heathland'])
Exclusion_area = Exclusion_area.rename_axis(index='descriptive names')

Inclusion_area = pd.DataFrame({"land classes" :[211, 212, 213, 231, 241, 242, 243, 321, 323, 324, 333],
                   "code value" : [12, 13, 14, 18, 19, 20, 21, 26, 28, 29, 32],}, 
                    index = ['Non-irrigated arable land', 'Permanently irrigated land', 'Rice fields', 
                  'Pastures', 'Annual crops associated with permanent crops', 'Complex cultivation patterns', 
                  'Land principally occupied by agriculture with significant areas of natural vegetation', 
                  'Natural grasslands', 'Sclerophyllous vegetation', 'Transitional woodland-shrub', 
                  'Sparsely vegetated areas'])
Inclusion_area = Inclusion_area.rename_axis(index='descriptive names')

# Defining function to plot the eligible area on the country map based on excluded & included areas:
def plot_area(masked, transform, shape):
    fig, ax = plt.subplots(figsize=(5, 5))
    ax = show(masked, transform=transform, cmap='Greens', vmin=0, ax=ax)
    shape.plot(ax=ax, edgecolor='k', color='None', linewidth=2)

# Defining a function to calculate the eligible area portion of the country in percentage based on excluded and included areas:
def eligible_area(masked, excluder, shape):
    return masked.sum() * float(excluder.res**2) / shape.geometry.area.sum() * 100

####################################################################################################

# Excluding the unsuitable area for onshore wind technology:

excluder_onshore_wind = ExclusionContainer()

excluder_onshore_wind.add_geometry("ne_10m_roads.gpkg", buffer=300)                 # 300m distance to major roads
excluder_onshore_wind.add_geometry("ne_10m_airports (1).gpkg", buffer=10000)        # 10km distance to airports
excluder_onshore_wind.add_raster('WDPA_Oct2022_Public_shp-SVK.tif', crs=3035)       # no natural protection areas
excluder_onshore_wind.add_raster('PROBAV_LC100_global_v3_0_1_2019_nrt_Discrete_Classification_map.tif', codes=[1,2,3,4,5,6,11,40,41], buffer=1000, crs=3035)  #1000m distance to built-up areas

onshore_clases = [12, 13, 14, 18, 19, 20, 21, 26, 28, 29, 32]         # only on suitable land cover classes
excluder_onshore_wind.add_raster('PROBAV_LC100_global_v3_0_1_2019_nrt_Discrete_Classification_map.tif', codes=onshore_clases, crs=3035)

masked, transform = shape_availability(shape_svk, excluder_onshore_wind)

# Calculating the eligible area for onshore wind technology in percentage (%):
onshore_eligible_area = eligible_area(masked, excluder_onshore_wind, shape_svk)
print()
print(f'Almost {round(onshore_eligible_area,3)} % of Slovakia is eligible for onshore wind technology development')
print()

# Plotting the eligible area for onshore wind technology:
plot_area(masked, transform, shape_svk)
plt.text(0.5, -0.4,f'Almost {round(onshore_eligible_area,3)} % of Slovakia is eligible for onshore wind technology development', ha='center', va='center', transform=plt.gca().transAxes)

plt.savefig('on-shore wind power eligibility.png')

plt.show()

Almost 5.432 % of Slovakia is eligible for onshore wind technology development

# Excluding the unsuitable area for solar power plants:

excluder_solar = ExclusionContainer()
Solar_classes = [1, 2, 3, 4, 5, 6, 11, 18, 26, 27, 30, 31, 41]
excluder_solar.add_raster('PROBAV_LC100_global_v3_0_1_2019_nrt_Discrete_Classification_map.tif', codes=Solar_classes, crs=3035)
excluder_solar = ExclusionContainer()
excluder_solar.add_raster('WDPA_Oct2022_Public_shp-SVK.tif', crs=3035)

masked, transform = shape_availability(shape_svk, excluder_solar)

# Calculating the eligible area for solar development in percentage (%):
eligible_area_solar = eligible_area(masked, excluder_solar,shape_svk)
print()
print(f'Almost {round(eligible_area_solar,2)} % of Slovakia land is eligible for solar')
print()

# Plotting the eligible area for solar development:
plot_area(masked, transform, shape_svk)
plt.text(0.5, -0.4,f'Almost {round(eligible_area_solar,2)} % of Slovakia land is eligible for solar implementation.', ha='center', va='center', transform=plt.gca().transAxes)

plt.savefig('Solar power eligibility.png')

plt.show()

Almost 62.69 % of Slovakia land is eligible for solar

# Define the file path for the data, being downloaded from the ERA5 dataset:
data_file_path = "era5-2017-SVK.nc"

# defining the bound limits to create a spatial extent:
minx, miny, maxx, maxy = svk_regions.total_bounds
buffer = 0.25              # A buffer of 0.25 degrees has been considered as requested.

# Creating the cutout for the year 2017:
cutout = atlite.Cutout(
    path=data_file_path,
    module="era5",
    x=slice(minx - buffer, maxx + buffer),
    y=slice(miny - buffer, maxy + buffer),
    time="2017",
)

# Calling the function cutout.prepare() & initiate the download:
cutout.prepare()

# the downloaded data is accessible here below :
cutout.data

INFO:atlite.data:Cutout already prepared.

<xarray.Dataset>
Dimensions:           (x: 25, y: 10, time: 8760)
Coordinates:
  * x                 (x) float64 16.75 17.0 17.25 17.5 ... 22.25 22.5 22.75
  * y                 (y) float64 47.5 47.75 48.0 48.25 ... 49.25 49.5 49.75
  * time              (time) datetime64[ns] 2017-01-01 ... 2017-12-31T23:00:00
    lon               (x) float64 dask.array<chunksize=(25,), meta=np.ndarray>
    lat               (y) float64 dask.array<chunksize=(10,), meta=np.ndarray>
Data variables: (12/13)
    height            (y, x) float32 dask.array<chunksize=(10, 25), meta=np.ndarray>
    wnd100m           (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    wnd_azimuth       (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    roughness         (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    influx_toa        (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    influx_direct     (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    ...                ...
    albedo            (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    solar_altitude    (time, y, x) float64 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    solar_azimuth     (time, y, x) float64 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    temperature       (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    soil temperature  (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
    runoff            (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
Attributes:
    module:             era5
    prepared_features:  ['temperature', 'influx', 'wind', 'runoff', 'height']
    chunksize_time:     100
    Conventions:        CF-1.6
    history:            2024-01-31 15:41:18 GMT by grib_to_netcdf-2.25.1: /op...

array([16.75, 17.  , 17.25, 17.5 , 17.75, 18.  , 18.25, 18.5 , 18.75, 19.  ,
       19.25, 19.5 , 19.75, 20.  , 20.25, 20.5 , 20.75, 21.  , 21.25, 21.5 ,
       21.75, 22.  , 22.25, 22.5 , 22.75])

array([47.5 , 47.75, 48.  , 48.25, 48.5 , 48.75, 49.  , 49.25, 49.5 , 49.75])

array(['2017-01-01T00:00:00.000000000', '2017-01-01T01:00:00.000000000',
       '2017-01-01T02:00:00.000000000', ..., '2017-12-31T21:00:00.000000000',
       '2017-12-31T22:00:00.000000000', '2017-12-31T23:00:00.000000000'],
      dtype='datetime64[ns]')

PandasIndex(Index([16.75,  17.0, 17.25,  17.5, 17.75,  18.0, 18.25,  18.5, 18.75,  19.0,
       19.25,  19.5, 19.75,  20.0, 20.25,  20.5, 20.75,  21.0, 21.25,  21.5,
       21.75,  22.0, 22.25,  22.5, 22.75],
      dtype='float64', name='x'))

PandasIndex(Index([47.5, 47.75, 48.0, 48.25, 48.5, 48.75, 49.0, 49.25, 49.5, 49.75], dtype='float64', name='y'))

PandasIndex(DatetimeIndex(['2017-01-01 00:00:00', '2017-01-01 01:00:00',
               '2017-01-01 02:00:00', '2017-01-01 03:00:00',
               '2017-01-01 04:00:00', '2017-01-01 05:00:00',
               '2017-01-01 06:00:00', '2017-01-01 07:00:00',
               '2017-01-01 08:00:00', '2017-01-01 09:00:00',
               ...
               '2017-12-31 14:00:00', '2017-12-31 15:00:00',
               '2017-12-31 16:00:00', '2017-12-31 17:00:00',
               '2017-12-31 18:00:00', '2017-12-31 19:00:00',
               '2017-12-31 20:00:00', '2017-12-31 21:00:00',
               '2017-12-31 22:00:00', '2017-12-31 23:00:00'],
              dtype='datetime64[ns]', name='time', length=8760, freq=None))

# Which data is contained in the cutout:
cutout

<Cutout "era5-2017-SVK">
 x = 16.75 ⟷ 22.75, dx = 0.25
 y = 47.50 ⟷ 49.75, dy = 0.25
 time = 2017-01-01 ⟷ 2017-12-31, dt = H
 module = era5
 prepared_features = ['height', 'wind', 'influx', 'temperature', 'runoff']

# Included weather variables list:

cutout.prepared_features

module  feature    
era5    height                   height
        wind                    wnd100m
        wind                wnd_azimuth
        wind                  roughness
        influx               influx_toa
        influx            influx_direct
        influx           influx_diffuse
        influx                   albedo
        influx           solar_altitude
        influx            solar_azimuth
        temperature         temperature
        temperature    soil temperature
        runoff                   runoff
dtype: object

# Creating availability matrixes for solar:
?cutout.availabilitymatrix
solar_AM = cutout.availabilitymatrix(shape_svk, excluder_solar)
#display(solar_AM)

# Creating availability matrixes for onshore wind:
onshore_AM = cutout.availabilitymatrix(shape_svk, excluder_onshore_wind)
#display(onshore_AM)

cap_per_sqkm = 3          # [MW/km2]  For both wind and solar, a deployment density of 3 MW/km2 has been assumed.

# Calculating the area :
country_area = cutout.grid.set_index(["y", "x"]).to_crs(3035).area / 1e6
country_area = xr.DataArray(country_area, dims=("spatial"))
#display(country_area)

Compute availability matrix: 100%|███████████████████████████████████████████████| 8/8 [00:01<00:00,  4.56 gridcells/s]
Compute availability matrix: 100%|███████████████████████████████████████████████| 8/8 [01:43<00:00, 12.97s/ gridcells]

Signature:
cutout.availabilitymatrix(
    shapes,
    excluder,
    nprocesses=None,
    disable_progressbar=False,
)
Docstring:
Compute the eligible share within cutout cells in the overlap with shapes.

For parallel calculation (nprocesses not None) the excluder must not be
initialized and all raster references must be strings. Otherwise processes
are colliding when reading from one common rasterio.DatasetReader.

Parameters
----------
cutout : atlite.Cutout
    Cutout which the availability matrix is aligned to.
shapes : geopandas.Series/geopandas.DataFrame
    Geometries for which the availabilities are calculated.
excluder : atlite.gis.ExclusionContainer
    Container of all meta data or objects which to exclude, i.e.
    rasters and geometries.
nprocesses : int, optional
    Number of processes to use for calculating the matrix. The paralle-
    lization can heavily boost the calculation speed. The default is None.
disable_progressbar: bool, optional
    Disable the progressbar if nprocesses is not None. Then the `map`
    function instead of the `imap` function is used for the multiprocessing
    pool. This speeds up the calculation.

Returns
-------
availabilities : xr.DataArray
    DataArray of shape (|shapes|, |y|, |x|) containing all the eligible
    share of cutout cell (x,y) in the overlap with shape i.

Notes
-----
The rasterio (or GDAL) average downsampling returns different results
dependent on how the target raster (the cutout raster) is spanned.
Either it is spanned from the top left going downwards,
e.g. Affine(0.25, 0, 0, 0, -0.25, 50), or starting in the
lower left corner and going up, e.g. Affine(0.25, 0, 0, 0, 0.25, 50).
Here we stick to the top down version which is why we use
`cutout.transform_r` and flipping the y-axis in the end.
File:      c:\users\mahtab\anaconda3\envs\esm-2023\lib\site-packages\atlite\gis.py
Type:      method

# Creating capacity matrix for onshore wind:
capacity_matrix_onshore = onshore_AM.stack(spatial=["y", "x"]) * country_area * cap_per_sqkm
#display(capacity_matrix_onshore)

# Wind profile:
wind_onshore = cutout.wind(
    atlite.windturbines.Vestas_V112_3MW,
    matrix=capacity_matrix_onshore,
    index= svk_regions.index,
    per_unit=True,
)

svk_wind_onshore_2017_potential_time_series = wind_onshore.to_pandas()
display(svk_wind_onshore_2017_potential_time_series.head(15))

INFO:atlite.convert:Convert and aggregate 'wind'.

# Plotting the on-shore wind potential time series for Jan-2017: 

plt.figure(figsize=(15, 6))  

region_colors = {
    'Banskobystrický': 'b',  
    'Bratislavský': 'g',  
    'Košický': 'r',  
    'Nitriansky': 'y',  
    'Prešovský': 'brown',  
    'Trenčiansky': 'pink',  
    'Trnavský': 'k', 
    'Žilinský': 'c', 
}

for region, color in region_colors.items():
    svk_wind_onshore_2017_potential_time_series.loc['2017-01', region].plot(color=color, label=region)

plt.ylabel('On-shore Wind Potential [p.u.]')  
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] of all regions in the January, 2017')
plt.savefig('On-shore Wind Potential-January.png')
plt.legend()
plt.grid(True)
plt.show()


################################################

# Plotting the on-shore wind potential time series: ( Banskobystrický)

plt.figure(figsize=(12, 6))  
svk_wind_onshore_2017_potential_time_series.iloc[:,0].plot(color='red')
plt.ylabel('On-shore Wind Potential [p.u.]')  
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] for Banskobystrický region in the country of Slovokia')
plt.savefig('On-shore Wind Potential-Banskobystrický.png')

plt.show()

##################################################
# Plotting the on-shore wind potential time series: ( Bratislavský)

plt.figure(figsize=(12, 6))  
svk_wind_onshore_2017_potential_time_series.iloc[:,1].plot(color='red')
plt.ylabel('On-shore Wind Potential [p.u.]')  
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] for Bratislavský region in the country of Slovokia')
plt.savefig('On-shore Wind Potential-Bratislavský.png')

plt.show()

################################################

# Plotting the on-shore wind potential time series: ( Košický)

plt.figure(figsize=(12, 6))  
svk_wind_onshore_2017_potential_time_series.iloc[:,2].plot(color='green')
plt.ylabel('On-shore Wind Potential [p.u.]')  
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] for Košický region in the country of Slovokia')
plt.savefig('On-shore Wind Potential-Košický.png')

plt.show()


################################################

# Plotting the on-shore wind potential time series: ( Nitriansky)

plt.figure(figsize=(12, 6))  
svk_wind_onshore_2017_potential_time_series.iloc[:,3].plot(color='grey')
plt.ylabel('On-shore Wind Potential [p.u.]')  
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] for Nitriansky region in the country of Slovokia')
plt.savefig('On-shore Wind Potential-Nitriansky.png')

plt.show()

# Creating capacity matrix for solar:
capacity_matrix_solar = solar_AM.stack(spatial=['y', 'x']) * country_area * cap_per_sqkm
#display(capacity_matrix_solar)

# Solar PV profile:
pv = cutout.pv(
    panel=atlite.solarpanels.CdTe,
    matrix=capacity_matrix_solar,
    orientation="latitude_optimal",
    index= svk_regions.index,
    per_unit=True,
)

svk_pv_2017_potential_time_series = pv.to_pandas()
svk_pv_2017_potential_time_series.tail(15)

INFO:atlite.convert:Convert and aggregate 'pv'.

[########################################] | 100% Completed | 5.36 s

# Plotting the solar potential time series for Jan-2017: 
plt.figure(figsize=(15, 6))  

for region, color in region_colors.items():
    pv.loc['2017-01', region].plot(color=color, label=region)

plt.ylabel('On-shore Wind Potential [p.u.]')  
plt.ylim(0, 1)
plt.title('Solar Potential [p.u.] of all regions in the January, 2017')
plt.savefig('Solar Potential-January.png')
plt.legend()
plt.grid(True)
plt.show()

##########################################################


# Plotting the solar potential time series: (only Banskobystrický)

plt.figure(figsize=(12, 6)) 
pv.to_pandas().iloc[:,0].plot(color='blue')
plt.ylabel('Solar Potential [p.u.]')  
plt.ylim(0, 1)

plt.title('Solar Potential [p.u.] for Banskobystrický region in the country of Slovokia')
plt.savefig('Solar Potential-Banskobystrický.png')

plt.show()

##################################################

# Plotting the solar potential time series: (only Bratislavský)

plt.figure(figsize=(12, 6)) 
pv.to_pandas().iloc[:,1].plot(color='red')
plt.ylabel('Solar Potential [p.u.]')  
plt.ylim(0, 1)

plt.title('Solar Potential [p.u.] for Bratislavský region in the country of Slovokia')
plt.savefig('Solar Potential-Bratislavský.png')

plt.show()


######################################################


# Plotting the solar potential time series: (only Košický)

plt.figure(figsize=(12, 6)) 
pv.to_pandas().iloc[:,2].plot(color='green')
plt.ylabel('Solar Potential [p.u.]')  
plt.ylim(0, 1)

plt.title('Solar Potential [p.u.] for Košický region in the country of Slovokia')
plt.savefig('Solar Potential-Košický.png')

plt.show()

################################################################

# Plotting the solar potential time series: (only Nitriansky)

plt.figure(figsize=(12, 6)) 
pv.to_pandas().iloc[:,3].plot(color='g')
plt.ylabel('Solar Potential [p.u.]')  
plt.ylim(0, 1)

plt.title('Solar Potential [p.u.] for Nitriansky region in the country of Slovokia')
plt.savefig('Solar Potential-Nitriansky.png')

plt.show()

from matplotlib.cm import ScalarMappable

## Building the Pypsa Model to minimize total annual system costs:

#display(gdf_svk.head(3))


n = pypsa.Network()       # defining a network

start = '2017-01-01 00:00:00'
end = '2017-12-31 23:00:00'

timestamps = pd.date_range(start=start, end=end, freq='H')

n.set_snapshots(timestamps)


# Defining the Buses:

regions_points = svk_regions.representative_point()

for index, row in svk_regions.iterrows():                    # A bus (Node) is defined for each region based on its geometry (latitude & Longitude)
    bus_name = index
    if bus_name in n.buses.index:
        n.remove("Bus", bus_name)
    region_point = regions_points.loc[index]  
    n.add("Bus", bus_name, x=region_point.x, y=region_point.y)

display(n.buses)

# !!!! FOR PRESENTATION: !!!
fig, ax = plt.subplots(figsize=(13, 9))
svk_regions.plot(ax=ax, facecolor='lightblue', edgecolor='black', linewidth=1.5)
regions_points.plot(ax=ax, color="black", markersize=20);
for bus_name, bus in n.buses.iterrows():
    ax.annotate(bus_name, (bus['x'], bus['y']), textcoords="offset points", xytext=(3,10), ha='center')
plt.title('regions of Slovokia as Nodes in the Model')
plt.savefig('Nodes.png')

plt.show()

##################################################################
# !!!! FOR PRESENTATION: !!!
fig, ax = plt.subplots(figsize=(13, 9))
svk_regions.plot(ax=ax, facecolor='lightblue', edgecolor='black', linewidth=1.5)
regions_points.plot(ax=ax, color="lightblue", markersize=0.5);


for bus_name, bus in n.buses.iterrows():
    ax.annotate(bus_name, (bus['x'], bus['y']), textcoords="offset points", xytext=(3,10), ha='center')

power_plants_scatter = gdf_svk.plot(
    ax=ax ,
    column = gdf_svk['primary_fuel'],
    markersize = gdf_svk['capacity_mw'] / 1e1 *4,
    legend= True,
    cmap='viridis',
)

plt.title('Power plants Capacities per region in Slovokia')
plt.savefig('Power plants capacities.png')

plt.show()

year = 2040    # A year for cost projections has been selected
url = f"https://raw.githubusercontent.com/PyPSA/technology-data/master/outputs/costs_{year}.csv"
tech_properties = pd.read_csv(url)
#display(tech_properties.head(10))

tech_properties_svk = tech_properties[
    (tech_properties['technology'] == 'nuclear') |
    (tech_properties['technology'] == 'CCGT') |
    (tech_properties['technology'] == 'coal') |
    (tech_properties['technology'] == 'hydro') |
    (tech_properties['technology'] == 'solar') |
    (tech_properties['technology'] == 'gas') |
    (tech_properties['technology'] == 'electrolysis') |
    (tech_properties['technology'] == 'fuel cell') |
    (tech_properties['technology'] == 'hydrogen storage underground') |
    (tech_properties['technology'] == 'battery inverter') |
    (tech_properties['technology'] == 'battery storage')
    
]
tech_properties_svk['technology'] = tech_properties_svk['technology'].replace('gas', 'CCGT')
display(tech_properties_svk)

# Dataframe of the technologies:

df_CCGT = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'CCGT']
df_CCGT.set_index('parameter', inplace=True)
#display(df_CCGT)

df_solar = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'solar']
df_solar.set_index('parameter', inplace=True)
df_solar.loc['fuel'] = ['solar', 0, None, None, None]
df_solar.loc['efficiency'] = ['solar', 1, None, None, None]
#display(df_solar)

df_nuclear = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'nuclear']
df_nuclear.set_index('parameter', inplace=True)
#display(df_nuclear)

df_hydro = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'hydro']
df_hydro.set_index('parameter', inplace=True)
df_hydro.loc['fuel'] = ['hydro', 0, None, None, None]
df_hydro.loc['VOM'] = ['hydro', 0, None, None, None]
#display(df_hydro)

df_coal = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'coal']
df_coal.set_index('parameter', inplace=True)
#display(df_coal)

df_electrolysis = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'electrolysis']
df_electrolysis.set_index('parameter', inplace=True)

df_fuel_cell = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'fuel cell']
df_fuel_cell.set_index('parameter', inplace=True)

df_hydrogen_storage_underground = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'hydrogen storage underground']
df_hydrogen_storage_underground.set_index('parameter', inplace=True)
display(df_electrolysis)

df_battery_inverter = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'battery inverter']
df_battery_inverter.set_index('parameter', inplace=True)

df_battery_storage = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'battery storage']
df_battery_storage.set_index('parameter', inplace=True)

# Calculating the marginal costs of the technologies:

marginalcost_Coal = df_coal.loc['fuel', 'value'] / df_coal.loc['efficiency', 'value'] + df_coal.loc['VOM', 'value']
marginalcost_Solar = df_solar.loc['fuel', 'value'] / df_solar.loc['efficiency', 'value'] + df_solar.loc['VOM', 'value']
marginalcost_Nuclear = df_nuclear.loc['fuel', 'value'] / df_nuclear.loc['efficiency', 'value'] + df_nuclear.loc['VOM', 'value']
marginalcost_Gas = df_CCGT.loc['fuel', 'value'] / df_CCGT.loc['efficiency', 'value'] + df_CCGT.loc['VOM', 'value']
marginalcost_Hydro = df_hydro.loc['fuel', 'value'] / df_hydro.loc['efficiency', 'value'] + df_hydro.loc['VOM', 'value']
marginalcost_Wind = 0


marginal_costs ={ 'Coal' : marginalcost_Coal,
                 'Solar' : marginalcost_Solar,
                 'Nuclear': marginalcost_Nuclear,
                 'Gas': marginalcost_Gas,
                 'Hydro' : marginalcost_Hydro,
                 'Wind' : marginalcost_Wind,
                }

# Calculating the capital costs of the technologies:

r = 0.07

annuity_coal = (r / (1 - (1 + r) ** -df_coal.loc['lifetime', 'value'])) * (df_coal.loc['investment', 'value']) * 1000
annuity_solar = (r / (1 - (1 + r) ** -df_solar.loc['lifetime', 'value'])) * (df_solar.loc['investment', 'value']) * 1000
annuity_hydro = (r / (1 - (1 + r) ** -df_hydro.loc['lifetime', 'value'])) * (df_hydro.loc['investment', 'value']) * 1000
annuity_nuclear = (r / (1 - (1 + r) ** -df_nuclear.loc['lifetime', 'value'])) * (df_nuclear.loc['investment', 'value']) * 1000
annuity_CCGT = (r / (1 - (1 + r) ** -df_CCGT.loc['lifetime', 'value'])) * (df_CCGT.loc['investment', 'value']) * 1000
annuity_electrolysis = (r / (1 - (1 + r) ** -df_electrolysis.loc['lifetime', 'value'])) * (df_electrolysis.loc['investment', 'value']) * 1000
annuity_fuel_cell = (r / (1 - (1 + r) ** -df_fuel_cell.loc['lifetime', 'value'])) * (df_fuel_cell.loc['investment', 'value']) * 1000
annuity_hydrogen_storage_underground = (r / (1 - (1 + r) ** -df_hydrogen_storage_underground.loc['lifetime', 'value'])) * (df_hydrogen_storage_underground.loc['investment', 'value']) * 1000
annuity_battery_inverter = (r / (1 - (1 + r) ** -df_battery_inverter.loc['lifetime', 'value'])) * (df_battery_inverter.loc['investment', 'value']) * 1000
annuity_battery_storage = (r / (1 - (1 + r) ** -df_battery_storage.loc['lifetime', 'value'])) * (df_battery_storage.loc['investment', 'value']) * 1000

FOM_coal = (df_coal.loc['FOM', 'value'] / 100) * df_coal.loc['investment', 'value'] *1000
FOM_solar = (df_solar.loc['FOM', 'value'] / 100) * df_solar.loc['investment', 'value'] *1000
FOM_nuclear = (df_nuclear.loc['FOM', 'value'] / 100) * df_nuclear.loc['investment', 'value'] *1000
FOM_hydro = (df_hydro.loc['FOM', 'value'] / 100) * df_hydro.loc['investment', 'value'] *1000
FOM_CCGT = (df_CCGT.loc['FOM', 'value'] / 100) * df_CCGT.loc['investment', 'value'] *1000
FOM_electrolysis = (df_electrolysis.loc['FOM', 'value'] / 100) * df_electrolysis.loc['investment', 'value'] *1000
FOM_fuel_cell = (df_fuel_cell.loc['FOM', 'value'] / 100) * df_fuel_cell.loc['investment', 'value'] *1000
FOM_battery_inverter = (df_battery_inverter.loc['FOM', 'value'] / 100) * df_battery_inverter.loc['investment', 'value'] *1000


capitalcost_coal = annuity_coal + FOM_coal
capitalcost_solar = annuity_solar + FOM_solar
capitalcost_nuclear = annuity_nuclear + FOM_nuclear
capitalcost_hydro = annuity_hydro + FOM_hydro
capitalcost_CCGT = annuity_CCGT + FOM_CCGT
capitalcost_electrolysis= annuity_electrolysis + FOM_electrolysis
capitalcost_fuel_cell= annuity_fuel_cell + FOM_fuel_cell
capitalcost_hydrogen_storage_underground = annuity_hydrogen_storage_underground 
capitalcost_battery_storage = annuity_battery_storage
capitalcost_battery_inverter= annuity_battery_inverter + FOM_battery_inverter

display(capitalcost_battery_storage)

7575.121930044452

## Defining the Generators:

power_plants_df_with_regions = gpd.sjoin(gdf_svk, svk_regions, how="left", op="within")  # spatial join to specify in which region each power plant is located
#display(power_plants_df_with_regions.head(5))

aggregated_power_plants = power_plants_df_with_regions.groupby(['primary_fuel', 'index_right']).first()  # grouping the power plant data based on technology and region

#display(aggregated_power_plants)

# Adding the Carriers & Corresponding emissions:
emissions = {
    'Coal': 0.34,
    'Gas': 0.2,
    'Hydro': 0,
    'Nuclear': 0,
    'Solar': 0,
    'Wind': 0
}

n.madd(
    "Carrier",
    ["Coal", "Gas", "Hydro", "Wind", "Solar", "Nuclear"],
    co2_emissions=emissions,
)


# Adding the generators to the network

for index, row in aggregated_power_plants.iterrows():
    generator_name = f"{index[0]}_{index[1]}"       # Specifying a name for each generator
    
    if generator_name in n.generators.index:        # Removing already created generators to avoid errors regarding duplication of generators name
        n.remove("Generator", generator_name)
        
    p_nom = row['capacity_mw']

    if index[0]=='Hydro':
        p_max_pu = (row['estimated_generation_gwh_2017'] * 1000) / (row['capacity_mw'] * 365 * 24 )
    else:
        p_max_pu = 1
        
    n.add("Generator",                           # Adding generators
          generator_name, 
          bus=index[1],
          p_nom_extendable = False,              # The Generators are not extendable
          capital_cost = 0,                      # The Generators have no capital costs
          p_nom=p_nom,
          marginal_cost = marginal_costs[index[0]],
          p_max_pu = p_max_pu,
          carrier=index[0],
         )
    
    for generator_name in n.generators.index:                         # Disregarding existing wind and solar capacities
        if "Wind" in generator_name or "Solar" in generator_name:
            n.remove("Generator", generator_name)
        
display(n.generators.carrier)

Generator
Coal_Košický                Coal
Coal_Trenčiansky            Coal
Gas_Bratislavský             Gas
Gas_Nitriansky               Gas
Gas_Trnavský                 Gas
Hydro_Banskobystrický      Hydro
Hydro_Bratislavský         Hydro
Hydro_Košický              Hydro
Hydro_Trenčiansky          Hydro
Hydro_Trnavský             Hydro
Hydro_Žilinský             Hydro
Nuclear_Nitriansky       Nuclear
Nuclear_Trnavský         Nuclear
Name: carrier, dtype: object

load_time_series = pd.read_csv('load.csv', index_col=0, parse_dates=True)   # Reading the load file
load_data_series_svk = load_time_series[['SK']].copy()                      # Cutting up the Slovakia part of the time series
#display(load_data_series_svk.head(10))

new_index = load_data_series_svk.index.strftime('%Y-%m-%d %H:%M:%S').str.replace('2013', '2017')
load_data_series_svk.index = pd.to_datetime(new_index)
load_data_series_svk = load_data_series_svk.rename_axis('Time')


#display(load_data_series_svk)

# Defining a dictionary including regions population [ref.: https://www.citypopulation.de/en/slovakia/cities/]:
cities_population = {'Banskobystrický':617777,
                     'Bratislavský' : 728370,
                     'Košický': 779505,
                     'Nitriansky': 670696,
                     'Prešovský': 808090,
                     'Trenčiansky': 570675,
                     'Trnavský': 565573,
                     'Žilinský': 688106}

total_population_skv = sum(cities_population.values())    # Calculating the total population of Slovokia

# Adding extra columns as region's load value time series to the primary load time series:
for region in cities_population.keys():
    load_data_series_svk ['skv_load_per_capita'] = load_data_series_svk['SK'] / total_population_skv    # calculating the load per capita
    load_data_series_svk [f'load_{region}'] = load_data_series_svk ['skv_load_per_capita'] * cities_population [region]   # Calculating the load on each region based on the region's population

#display(load_data_series_svk)


for index , bus in n.buses.iterrows():
    Load_name = f'Demand_load_{index}'
    if (n.loads.index == Load_name).any():
        n.remove("Load", Load_name)
    p_set_data = load_data_series_svk[f'load_{index}'] 

    n.add("Load",                           # Adding Electricity loads for each bus (region)
        Load_name, 
        bus=index,
        p_set=p_set_data,             
        carrier="electricity",                     
        )

display(n.loads_t.p_set)

ax = n.loads_t.p_set.iloc[:24].plot(figsize=(10, 6),ylim=[0, 800], ylabel="Load [MW]", xlabel="Time", title='Electricity Demand Load during a single day in different regions of Slovokia (01.01.2017)')
plt.savefig('electricity_demand_load.png')
plt.show()
#n.snapshots

# Defining regions area as a dictionary[km^2]:
Regions_area = {
                'Banskobystrický' :  9455,
                'Bratislavský' : 2053,
                'Košický': 6753,
                'Nitriansky': 6343,
                'Prešovský' : 8974 ,
                'Trenčiansky': 4502 ,
                'Trnavský': 4145,
                'Žilinský': 6809,
}

# Adding solar and on-shore wind generators per region:
n.buses

for index , bus in n.buses.iterrows():
    solar_generator_name = f'Solar_{index}'
    if solar_generator_name in n.generators.index:
        n.remove("Generator", solar_generator_name)   
    n.add("Generator",                           # Adding solar generators
        solar_generator_name, 
        bus=index,
        p_nom_extendable = False,              # The Generators are not extendable
        capital_cost = 0,                      # The Generators have no capital costs
        marginal_cost = marginalcost_Solar,
        p_nom= Regions_area[index] * cap_per_sqkm,
        p_max_pu = svk_pv_2017_potential_time_series[index],
        carrier='Solar')

    wind_generator_name = f'Wind_{index}'

    if wind_generator_name in n.generators.index:
        n.remove("Generator", wind_generator_name)
    n.add("Generator",                           # Adding wind generators
        wind_generator_name, 
        bus=index,
        p_nom_extendable = False,              # The Generators are not extendable
        capital_cost = 0,                      # The Generators have no capital costs
        marginal_cost = marginalcost_Wind,
        p_nom= Regions_area[index] * cap_per_sqkm,
        p_max_pu = svk_wind_onshore_2017_potential_time_series[index],
        carrier='Wind',
         )


display( n.generators.carrier.map(n.carriers.co2_emissions))
display(n.generators_t.p_max_pu)

Generator
Coal_Košický             0.34
Coal_Trenčiansky         0.34
Gas_Bratislavský         0.20
Gas_Nitriansky           0.20
Gas_Trnavský             0.20
Hydro_Banskobystrický    0.00
Hydro_Bratislavský       0.00
Hydro_Košický            0.00
Hydro_Trenčiansky        0.00
Hydro_Trnavský           0.00
Hydro_Žilinský           0.00
Nuclear_Nitriansky       0.00
Nuclear_Trnavský         0.00
Solar_Banskobystrický    0.00
Wind_Banskobystrický     0.00
Solar_Bratislavský       0.00
Wind_Bratislavský        0.00
Solar_Košický            0.00
Wind_Košický             0.00
Solar_Nitriansky         0.00
Wind_Nitriansky          0.00
Solar_Prešovský          0.00
Wind_Prešovský           0.00
Solar_Trenčiansky        0.00
Wind_Trenčiansky         0.00
Solar_Trnavský           0.00
Wind_Trnavský            0.00
Solar_Žilinský           0.00
Wind_Žilinský            0.00
Name: carrier, dtype: float64

# Comparing solar and wind potential in January 2017:

n.generators_t.p_max_pu.groupby(n.generators.carrier, axis=1).mean().iloc[:744].plot(ylabel="p.u.", figsize=(15, 6), color=['y', 'b'])
plt.savefig('Mean renewable energy capacity of implemented Generators-Jan-2017.png')
plt.show()

# Comparing solar and wind potential in 1.January 2017 (as a single day):
n.generators_t.p_max_pu.groupby(n.generators.carrier, axis=1).mean().iloc[:24].plot(ylabel="p.u.", figsize=(15, 6), color=['y', 'b'])
plt.savefig('Mean renewable energy capacity of implemented Generators-one day-1.Jan-2017.png')
plt.show()

# Difining a list of tuples, including neighboring regions to create transmission lines between them:
neighbor_buses = [('Bratislavský','Trnavský'),
                  ('Trnavský','Nitriansky'),
                  ('Trnavský','Trenčiansky'),
                  ('Nitriansky','Trenčiansky'),
                  ('Nitriansky','Banskobystrický'),
                  ('Banskobystrický','Trenčiansky'),
                  ('Žilinský','Trenčiansky'),
                  ('Banskobystrický','Žilinský'),
                  ('Prešovský','Žilinský'),
                  ('Prešovský','Banskobystrický'),
                  ('Košický','Banskobystrický'),
                  ('Košický','Prešovský')]

for region1, region2 in neighbor_buses:  
    
    distance = np.sqrt((regions_points[region1].x - regions_points[region2].x) ** 2 +    # Calculating the distance for each transmission line based on the representative points of buses
                       (regions_points[region1].y - regions_points[region2].y) ** 2)   
    length = 1.5 * distance   # As requested, the length of 1.5 times the crow-fly distance between the regions’ representative points
    cost = 400 * length       # Assumption for costs of 400€/MW/km

    # Adding bidirectional links between neighboring regions
    TL1_name = f"TL1_{region1}_{region2}"
    if TL1_name in n.links.index:
        n.remove("Link", TL1_name)
    n.add("Link",
          TL1_name,
          bus0=region1,
          bus1=region2,
          capital_cost=cost,
          p_max_pu=1.0,         # no transmission losses
          efficiency=1.0,       # no transmission losses
          length=length,
         )

    TL2_name = f"TL2_{region2}_{region1}"
    if TL2_name in n.links.index:
        n.remove("Link", TL2_name)
    n.add("Link",
          TL2_name,
          bus0=region2,
          bus1=region1,
          capital_cost=cost,
          p_max_pu=1.0,           # no transmission losses
          efficiency=1.0,         # no transmission losses
          length=length,
         )

display(n.links)

# Adding Battery storage units:

for bus in n.buses.index:
    battery_storage_U_name = f"Battery_storage_{bus}"
    if battery_storage_U_name in n.storage_units.index:
        n.remove("StorageUnit", battery_storage_U_name)
    n.add(
        "StorageUnit",
        battery_storage_U_name,
        max_hours= 6,
        capital_cost = 6 * capitalcost_battery_storage + capitalcost_battery_inverter,
        bus = bus,
        marginal_cost =0,
        carrier = "battery storage",
        efficiency_dispatch = df_battery_inverter.loc['efficiency', 'value'],
        efficiency_store = df_battery_inverter.loc['efficiency', 'value'],
        cyclic_state_of_charge = True,
        p_nom_extendable = True,
  
    )

n.storage_units

# Adding hydrogen storage units

for bus in n.buses.index:
    storage_U_name = f"hydrogen_storage_{bus}"
    if storage_U_name in n.storage_units.index:
        n.remove("StorageUnit", storage_U_name)
    n.add(
        "StorageUnit",
        storage_U_name,
        max_hours= 336,
        capital_cost = 336 * capitalcost_hydrogen_storage_underground+capitalcost_fuel_cell+capitalcost_electrolysis,
        bus = bus,
        marginal_cost = 0,
        carrier = "hydrogen storage underground",
        efficiency_dispatch = df_fuel_cell.loc['efficiency', 'value'],
        efficiency_store = df_electrolysis.loc['efficiency', 'value'],
        cyclic_state_of_charge = True,
        p_nom_extendable = True,   
    )

display(n.storage_units)

# Adding shedding generator:

n.madd(
            "Generator",
            n.buses.index,
            suffix=" load shedding",
            bus=n.buses.index,
            carrier="load",
            sign=1e-3,  # Adjust sign to measure p and p_nom in kW instead of MW
            marginal_cost=1000,  # Eur/kWh
            p_nom=1e9,  # kW
)

Index(['Banskobystrický load shedding', 'Bratislavský load shedding',
       'Košický load shedding', 'Nitriansky load shedding',
       'Prešovský load shedding', 'Trenčiansky load shedding',
       'Trnavský load shedding', 'Žilinský load shedding'],
      dtype='object', name='Bus')

# Calculating the CO2 emissions for each region based on the resulted power output:

e = (
    n.generators_t.p
    / n.generators.efficiency
    * n.generators.carrier.map(n.carriers.co2_emissions)
)
e.tail(20)

# Calculating the CO2 emissions in total (for the whole country)[tonnes/MWh thermal]:
display(e.sum().sum())

# Adding global constraint for CO2 emission 
n.add(
    "GlobalConstraint",
    "emission_limit",
    carrier_attribute="co2_emissions",
    sense="<=",
    constant=e.sum().sum() * 1,               # no reduction constraints in CO2 emission 
)


# Solving the Model:

n.optimize(solver_name="highs")

n.generators_t.p    # Retrieveing the power output

0.0

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:13<00:00,  1.08it/s]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:02<00:00,  2.52it/s]
INFO:linopy.io: Writing time: 17.16s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 4.59e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.

# Regional distribution of loads in Slovakia


load = n.loads_t.p_set.sum(axis=0).groupby(n.loads.bus).sum()
display(load)

fig = plt.figure(figsize=(15, 6))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.EqualEarth())

# Plotting the network with bus sizes proportional to the load
n.plot(
    ax=ax,
    bus_sizes=load / 2e9,  
    geomap=True, 
    color_geomap=True,  
    projection=ccrs.EqualEarth()
)
svk_regions.plot(ax=ax, facecolor='lightblue', edgecolor='black', linewidth=1.5, transform=ccrs.PlateCarree())
ax.set_extent([minx, maxx, miny, maxy], crs=ccrs.PlateCarree())  
plt.title('Regional Distribution of Load')
plt.show()

bus
Banskobystrický    3.480420e+06
Bratislavský       4.103476e+06
Košický            4.391560e+06
Nitriansky         3.778554e+06
Prešovský          4.552601e+06
Trenčiansky        3.215057e+06
Trnavský           3.186314e+06
Žilinský           3.876638e+06
dtype: float64

carrier_to_exclude = 'load'                # Excluding the "load" carriers, which are created during the Shedding Generator defining
load_excluded_generators = n.generators[n.generators.carrier != carrier_to_exclude]

dispatch_by_carrier = n.generators_t.p.groupby(load_excluded_generators.carrier, axis=1).sum().div(1e3)
#display(p_by_carrier.head(20))

fig, ax = plt.subplots(figsize=(15, 8))

dispatch_by_carrier.iloc[:720].plot(
    kind="area",
    ax=ax,
    linewidth=0,
    cmap="tab20b",
)

ax.legend(ncol=5, loc="upper left", frameon=False)

ax.set_ylabel("GW")
ax.set_title('Hourly dispatch per technology- JAN 2017')
ax.set_ylim(0, 9);

#########################################################
fig, ax = plt.subplots(figsize=(15, 8))

dispatch_by_carrier.iloc[:24].plot(
    kind="area",
    ax=ax,
    linewidth=0,
    cmap="tab20b",
)

ax.legend(ncol=5, loc="upper left", frameon=False)

ax.set_ylabel("GW")
ax.set_title('Hourly dispatch per technology through a single day- 1.JAN 2017')

ax.set_ylim(0, 9);

# The aggregate dispatch of the pumped hydro storage units and the state of charge throughout the day:

fig, ax = plt.subplots(figsize=(15, 8))

battery_storage_units = n.storage_units[n.storage_units.carrier == 'battery storage']

p_battery_storage = n.storage_units_t.p[battery_storage_units.index].sum(axis=1).div(1e3)
state_of_charge_battery = n.storage_units_t.state_of_charge[battery_storage_units.index].sum(axis=1).div(1e3)

hydro_storage_units = n.storage_units[n.storage_units.carrier == 'hydrogen storage underground']
p_hydro_storage = n.storage_units_t.p[hydro_storage_units.index].sum(axis=1).div(1e3)
state_of_charge_hydro = n.storage_units_t.state_of_charge[hydro_storage_units.index].sum(axis=1).div(1e3)

total_storage_region1 = p_hydro_storage + p_battery_storage
total_SOC_region1 = state_of_charge_hydro +  state_of_charge_battery
#display(p_hydro_storage.iloc[:24])
total_storage_region1.iloc[:24].plot(label=" dispatch [GW]", ax=ax)
total_SOC_region1.iloc[:24].plot(label=" state of charge [GWh]", ax=ax)

ax.grid()
ax.legend()
ax.set_ylabel("MWh or MW")
ax.set_title('Total Storage: Dispatch and State of Charge throughout the day in Banskobystrický ')

plt.show()
###############################################################

# The aggregate dispatch of the pumped hydro storage units and the state of charge throughout a Month (January):

fig, ax = plt.subplots(figsize=(15, 8))


total_storage_region1.iloc[:720].plot(label=" dispatch [GW]", ax=ax)
total_SOC_region1.iloc[:720].plot(label=" state of charge [GWh]", ax=ax)

ax.grid()
ax.legend()
ax.set_ylabel("MWh or MW")
ax.set_title('Total Storage: Dispatch and State of Charge throughout a month in Banskobystrický ')

plt.show()



##############################################################
# the aggregate dispatch of the pumped hydro storage units and the state of charge throughout the whole year:

fig, ax = plt.subplots(figsize=(15, 8))


total_storage_region1.plot(label=" dispatch [GW]", ax=ax)
total_SOC_region1.plot(label=" state of charge [GWh]", ax=ax)

ax.grid()
ax.legend()
ax.set_ylabel("MWh or MW")
ax.set_title('Total Storage: Dispatch and State of Charge throughout the whole year in Banskobystrický ')

plt.show()

# Locational marginal prices:


fig = plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.EqualEarth())

norm = plt.Normalize(vmin=0, vmax=50)  # €/MWh

n.plot(
    ax=ax,
    bus_colors=n.buses_t.marginal_price.mean(),
    bus_cmap="plasma",
    bus_norm=norm,

)

cax = fig.add_axes([0.95, 0.15, 0.02, 0.7])  
cb = plt.colorbar(
    plt.cm.ScalarMappable(cmap="plasma", norm=norm),
    label="LMP [€/MWh]",
    shrink=0.6,
    cax=cax,  
)
ax.set_extent([minx, maxx, miny, maxy], crs=ccrs.PlateCarree())  

plt.show()

display(n.generators['marginal_cost'])
generation_cost = (n.generators_t.p * n.generators['marginal_cost']).sum(axis=1)

sorted_costs = generation_cost[np.argsort(-generation_cost)]    # Sorting generation costs in ascending order
cumulative_probability = np.arange(1, len(sorted_costs) + 1) / len(sorted_costs)   # Calculating the cumulative probability

# Plotting:
plt.figure(figsize=(10, 6))
plt.plot(cumulative_probability * 100, sorted_costs, linestyle='-')
plt.xlabel('Cumulative Probability (%)')
plt.ylabel(' price (€)')
plt.title('Cost Duration Curve')
plt.grid(True)
plt.show()

Generator
Coal_Košický                       31.431739
Coal_Trenčiansky                   31.431739
Gas_Bratislavský                   44.520508
Gas_Nitriansky                     44.520508
Gas_Trnavský                       44.520508
Hydro_Banskobystrický               0.000000
Hydro_Bratislavský                  0.000000
Hydro_Košický                       0.000000
Hydro_Trenčiansky                   0.000000
Hydro_Trnavský                      0.000000
Hydro_Žilinský                      0.000000
Nuclear_Nitriansky                 13.778923
Nuclear_Trnavský                   13.778923
Solar_Banskobystrický               0.010000
Wind_Banskobystrický                0.000000
Solar_Bratislavský                  0.010000
Wind_Bratislavský                   0.000000
Solar_Košický                       0.010000
Wind_Košický                        0.000000
Solar_Nitriansky                    0.010000
Wind_Nitriansky                     0.000000
Solar_Prešovský                     0.010000
Wind_Prešovský                      0.000000
Solar_Trenčiansky                   0.010000
Wind_Trenčiansky                    0.000000
Solar_Trnavský                      0.010000
Wind_Trnavský                       0.000000
Solar_Žilinský                      0.010000
Wind_Žilinský                       0.000000
Banskobystrický load shedding    1000.000000
Bratislavský load shedding       1000.000000
Košický load shedding            1000.000000
Nitriansky load shedding         1000.000000
Prešovský load shedding          1000.000000
Trenčiansky load shedding        1000.000000
Trnavský load shedding           1000.000000
Žilinský load shedding           1000.000000
Name: marginal_cost, dtype: float64

n.global_constraints.mu

GlobalConstraint
emission_limit   -28014.414082
Name: mu, dtype: float64

#Setting a constraint to  solver:
n.global_constraints.loc["emission_limit", "constant"] = 0

 n.optimize(solver_name="highs")

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:15<00:00,  1.06s/it]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:02<00:00,  2.76it/s]
INFO:linopy.io: Writing time: 18.9s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 4.59e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.

('ok', 'optimal')

n.generators_t.p

carrier_to_exclude = 'load'                # Excluding the "load" carriers, which are created during the Shedding Generator defining
load_excluded_generators = n.generators[n.generators.carrier != carrier_to_exclude]

dispatch_by_carrier = n.generators_t.p.groupby(load_excluded_generators.carrier, axis=1).sum().div(1e3)
#display(p_by_carrier.head(20))

fig, ax = plt.subplots(figsize=(15, 8))

dispatch_by_carrier.iloc[:720].plot(
    kind="area",
    ax=ax,
    linewidth=0,
    cmap="tab20b",
)

ax.legend(ncol=5, loc="upper left", frameon=False)

ax.set_ylabel("GW")
ax.set_title('Hourly dispatch per technology- JAN 2017')
ax.set_ylim(0, 9);

#########################################################
fig, ax = plt.subplots(figsize=(15, 8))

dispatch_by_carrier.iloc[:24].plot(
    kind="area",
    ax=ax,
    linewidth=0,
    cmap="tab20b",
)

ax.legend(ncol=5, loc="upper left", frameon=False)

ax.set_ylabel("GW")
ax.set_title('Hourly dispatch per technology through a single day- 1.JAN 2017')

ax.set_ylim(0, 9);

generation_cost = (n.generators_t.p * n.generators['marginal_cost']).sum(axis=1)

sorted_costs = generation_cost[np.argsort(-generation_cost)]   # Sorting generation costs in ascending order

cumulative_probability = np.arange(1, len(sorted_costs) + 1) / len(sorted_costs)   # Calculating the cumulative probability


# Ploting:
plt.figure(figsize=(10, 6))
plt.plot(cumulative_probability * 100, sorted_costs, linestyle='-')
plt.xlabel('Cumulative Probability (%)')
plt.ylabel(' price (€)')
plt.title('Cost Duration Curve')
plt.grid(True)
plt.show()

n.global_constraints.mu

GlobalConstraint
emission_limit   -28014.414082
Name: mu, dtype: float64

# Locational marginal prices:

fig = plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.EqualEarth())

norm = plt.Normalize(vmin=0, vmax=50)  # €/MWh

n.plot(
    ax=ax,
    bus_colors=n.buses_t.marginal_price.mean(),
    bus_cmap="plasma",
    bus_norm=norm,

)

cax = fig.add_axes([0.95, 0.15, 0.02, 0.7])  
cb = plt.colorbar(
    plt.cm.ScalarMappable(cmap="plasma", norm=norm),
    label="LMP [€/MWh]",
    shrink=0.6,
    cax=cax,  
)
ax.set_extent([minx, maxx, miny, maxy], crs=ccrs.PlateCarree())  

plt.show()

# Reduction wind and solar potential in 25% steps:

for reduction_step in np.arange(0.75, -0.01, -0.25):
    
    previou_step = reduction_step + 0.25
    reduction_factor =  reduction_step / previou_step
    
    for generator_name, generator_data in n.generators.iterrows():
        if generator_data['carrier'] in ['Solar', 'Wind']:
           n.generators.at[generator_name, 'p_nom'] *= reduction_factor
        
    n.optimize(solver_name="highs")
    print(f'{(1- round(reduction_step * 100))} % Reduction of Renewable Energy installable potential')
    n.optimize(solver_name="highs")
    print()
    print(f'CO2 Shadow price in case of {(1- round(reduction_step * 100))} % Reduction of Renewable Energy installable potential:')
    print()
    display(n.global_constraints.mu)
    print()
    ###################################################################################################
    # Dispatch power per technology :
    carrier_to_exclude = 'load'                # Excluding the "load" carriers, which are created during the Shedding Generator defining
    load_excluded_generators = n.generators[n.generators.carrier != carrier_to_exclude]
    
    dispatch_by_carrier = n.generators_t.p.groupby(load_excluded_generators.carrier, axis=1).sum().div(1e3)
    
    fig, ax = plt.subplots(figsize=(15, 8))
    
    dispatch_by_carrier.iloc[:720].plot(
        kind="area",
        ax=ax,
        linewidth=0,
        cmap="tab20b",
    )
    
    ax.legend(ncol=5, loc="upper left", frameon=False)
    
    ax.set_ylabel("GW")
    ax.set_title(f'Hourly dispatch per technology ({(1- round(reduction_step * 100))} % Reduction of RE potential)')
    ax.set_ylim(0, 9);
    
    #-----------------------------------------------------------
    fig, ax = plt.subplots(figsize=(15, 8))
    
    dispatch_by_carrier.iloc[:24].plot(
        kind="area",
        ax=ax,
        linewidth=0,
        cmap="tab20b",
    )
    
    ax.legend(ncol=5, loc="upper left", frameon=False)
    
    ax.set_ylabel("GW")
    ax.set_title(f'Hourly dispatch per technology through a single day ({(1- round(reduction_step * 100))} % Reduction of RE potential)')
    
    ax.set_ylim(0, 9);

    #######################################################################################
       
    generation_cost = (n.generators_t.p * n.generators['marginal_cost']).sum(axis=1)
    sorted_costs = generation_cost[np.argsort(-generation_cost)]

    cumulative_probability = np.arange(1, len(sorted_costs) + 1) / len(sorted_costs)
    

    plt.figure(figsize=(10, 6))
    plt.plot(cumulative_probability * 100, sorted_costs, linestyle='-')
    plt.xlabel('Cumulative Probability (%)')
    plt.ylabel(' price (€)')
    plt.title(f'Cost Duration Curve ({(1- round(reduction_step * 100))} % Reduction of RE potential)')
    plt.grid(True)
    plt.show()
    ##############################################################################
    fig = plt.figure(figsize=(7, 7))
    ax = plt.axes(projection=ccrs.EqualEarth())
    
    norm = plt.Normalize(vmin=0, vmax=50)  # €/MWh
    
    n.plot(
        ax=ax,
        bus_colors=n.buses_t.marginal_price.mean(),
        bus_cmap="plasma",
        bus_norm=norm,
    
    )
    
    cax = fig.add_axes([0.95, 0.15, 0.02, 0.7])  
    cb = plt.colorbar(
        plt.cm.ScalarMappable(cmap="plasma", norm=norm),
        label="LMP [€/MWh]",
        shrink=0.6,
        cax=cax,  
    )
    ax.set_extent([minx, maxx, miny, maxy], crs=ccrs.PlateCarree())  
    ax.set_title(f'Electricity price per region ({(1- round(reduction_step * 100))} % Reduction of RE potential)')

    plt.show()

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:13<00:00,  1.10it/s]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:02<00:00,  2.43it/s]
INFO:linopy.io: Writing time: 17.02s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 4.84e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.

-74 % Reduction of Renewable Energy installable potential

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:19<00:00,  1.27s/it]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:03<00:00,  1.77it/s]
INFO:linopy.io: Writing time: 23.68s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 4.84e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.

CO2 Shadow price in case of -74 % Reduction of Renewable Energy installable potential:

GlobalConstraint
emission_limit   -28016.106241
Name: mu, dtype: float64

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:41<00:00,  2.74s/it]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:08<00:00,  1.17s/it]
INFO:linopy.io: Writing time: 50.74s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 5.32e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.

-49 % Reduction of Renewable Energy installable potential

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:39<00:00,  2.64s/it]
Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:07<00:00,  1.08s/it]
INFO:linopy.io: Writing time: 48.54s
INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log.
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 954856 primals, 2119937 duals
Objective: 5.32e+08
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.

CO2 Shadow price in case of -49 % Reduction of Renewable Energy installable potential:

	country	country_long	name	gppd_idnr	capacity_mw	latitude	longitude	primary_fuel	other_fuel1	other_fuel2	...	estimated_generation_gwh_2013	estimated_generation_gwh_2014	estimated_generation_gwh_2015	estimated_generation_gwh_2016	estimated_generation_gwh_2017	estimated_generation_note_2013	estimated_generation_note_2014	estimated_generation_note_2015	estimated_generation_note_2016	estimated_generation_note_2017
19861	SVK	Slovakia	Badín	WKS0063339	6.0	48.2650	19.0560	Solar	NaN	NaN	...	8.29	8.91	8.90	8.58	8.36	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE	SOLAR-V1-NO-AGE
19862	SVK	Slovakia	Bohunice Nuclear Power Plant Slovakia	GEODB0002385	880.0	48.4960	17.6816	Nuclear	NaN	NaN	...	NaN	NaN	NaN	NaN	6805.78	NO-ESTIMATION	NO-ESTIMATION	NO-ESTIMATION	NO-ESTIMATION	CAPACITY-FACTOR-V1
19863	SVK	Slovakia	Cierny Vah Pumped Storage Hydroelectric Power ...	GEODB0042635	735.0	49.0088	19.9122	Hydro	NaN	NaN	...	2597.78	2210.18	2927.73	2216.24	2210.18	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1
19864	SVK	Slovakia	Cunovo Hydroelectric Power Plant Slovakia	GEODB0042641	24.0	48.0299	17.2254	Hydro	NaN	NaN	...	118.03	117.69	129.12	122.91	130.73	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1
19865	SVK	Slovakia	Dubnica nad Vahom Hydroelectric Power Plant Sl...	GEODB0042662	16.5	48.9638	18.1454	Hydro	NaN	NaN	...	73.11	55.75	63.33	52.24	55.80	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1	HYDRO-V1

NAME_1	Banskobystrický	Bratislavský	Košický	Nitriansky	Prešovský	Trenčiansky	Trnavský	Žilinský
time
2017-01-01 00:00:00	0.000003	0.000060	0.007087	0.000000e+00	0.030756	0.000335	0.000001	0.013306
2017-01-01 01:00:00	0.000000	0.000005	0.007612	0.000000e+00	0.034854	0.000475	0.000040	0.014443
2017-01-01 02:00:00	0.000000	0.000066	0.007142	1.130307e-06	0.039974	0.000691	0.000331	0.015357
2017-01-01 03:00:00	0.000000	0.000171	0.003622	5.924064e-07	0.041654	0.001257	0.000336	0.017744
2017-01-01 04:00:00	0.000000	0.000144	0.003500	0.000000e+00	0.047368	0.001554	0.000341	0.020749
2017-01-01 05:00:00	0.000000	0.000156	0.003398	0.000000e+00	0.054182	0.002014	0.000318	0.027655
2017-01-01 06:00:00	0.000000	0.000167	0.003533	4.355234e-08	0.062805	0.002464	0.000342	0.035510
2017-01-01 07:00:00	0.000000	0.000406	0.004653	6.033786e-07	0.067074	0.002595	0.000470	0.042285
2017-01-01 08:00:00	0.000000	0.000687	0.004192	2.199573e-06	0.066282	0.002300	0.001231	0.044154
2017-01-01 09:00:00	0.000000	0.000343	0.005305	1.851060e-07	0.058913	0.001037	0.000373	0.044615
2017-01-01 10:00:00	0.000019	0.000218	0.005946	1.088100e-06	0.053568	0.000092	0.000304	0.027705
2017-01-01 11:00:00	0.000000	0.000181	0.007191	1.481370e-06	0.049395	0.000112	0.000274	0.025161
2017-01-01 12:00:00	0.000000	0.000299	0.008346	1.112223e-04	0.067957	0.000216	0.000336	0.031990
2017-01-01 13:00:00	0.000000	0.000523	0.008210	1.796761e-04	0.085420	0.000869	0.000486	0.044031
2017-01-01 14:00:00	0.000106	0.001866	0.008984	3.887791e-03	0.107708	0.005368	0.007772	0.060552

	technology	parameter	value	unit	source	further description
29	CCGT	FOM	3.3006	%/year	Danish Energy Agency, technology_data_for_el_a...	05 Gas turb. CC, steam extract.: Fixed O&M
30	CCGT	VOM	4.1000	EUR/MWh	Danish Energy Agency, technology_data_for_el_a...	05 Gas turb. CC, steam extract.: Variable O&M
31	CCGT	c_b	2.1000	50oC/100oC	Danish Energy Agency, technology_data_for_el_a...	05 Gas turb. CC, steam extract.: Cb coefficient
32	CCGT	c_v	0.1500	50oC/100oC	Danish Energy Agency, technology_data_for_el_a...	05 Gas turb. CC, steam extract.: Cv coefficient
33	CCGT	efficiency	0.5900	per unit	Danish Energy Agency, technology_data_for_el_a...	05 Gas turb. CC, steam extract.: Electricity ...
34	CCGT	investment	815.0000	EUR/kW	Danish Energy Agency, technology_data_for_el_a...	05 Gas turb. CC, steam extract.: Nominal inve...
35	CCGT	lifetime	25.0000	years	Danish Energy Agency, technology_data_for_el_a...	05 Gas turb. CC, steam extract.: Technical li...
412	battery inverter	FOM	0.5400	%/year	Danish Energy Agency, technology_data_catalogu...	: Fixed O&M
413	battery inverter	efficiency	0.9600	per unit	Danish Energy Agency, technology_data_catalogu...	: Round trip efficiency DC
414	battery inverter	investment	100.0000	EUR/kW	Danish Energy Agency, technology_data_catalogu...	: Output capacity expansion cost investment
415	battery inverter	lifetime	10.0000	years	Danish Energy Agency, technology_data_catalogu...	: Technical lifetime
416	battery storage	investment	94.0000	EUR/kWh	Danish Energy Agency, technology_data_catalogu...	: Energy storage expansion cost investment
417	battery storage	lifetime	30.0000	years	Danish Energy Agency, technology_data_catalogu...	: Technical lifetime
579	coal	CO2 intensity	0.3361	tCO2/MWh_th	Entwicklung der spezifischen Kohlendioxid-Emis...	NaN
580	coal	FOM	1.3100	%/year	Lazard's levelized cost of energy analysis - v...	Calculated based on average of listed range, i...
581	coal	VOM	3.3278	EUR/MWh_e	Lazard's levelized cost of energy analysis - v...	Calculated based on average of listed range, i...
582	coal	efficiency	0.3300	p.u.	Lazard's levelized cost of energy analysis - v...	Calculated based on average of listed range, i...
583	coal	fuel	9.2743	EUR/MWh_th	DIW (2013): Current and propsective costs of e...	Based on IEA 2011 data, 99 USD/t.
584	coal	investment	3905.3074	EUR/kW_e	Lazard's levelized cost of energy analysis - v...	Higher costs include coal plants with CCS, but...
585	coal	lifetime	40.0000	years	Lazard's levelized cost of energy analysis - v...	NaN
707	electrolysis	FOM	4.0000	%/year	Danish Energy Agency, data_sheets_for_renewabl...	86 AEC 100 MW: Fixed O&M
708	electrolysis	efficiency	0.6532	per unit	Danish Energy Agency, data_sheets_for_renewabl...	86 AEC 100 MW: Hydrogen Output
709	electrolysis	efficiency-heat	0.1849	per unit	Danish Energy Agency, data_sheets_for_renewabl...	86 AEC 100 MW: - hereof recoverable for dist...
710	electrolysis	investment	384.9356	EUR/kW_e	Danish Energy Agency, data_sheets_for_renewabl...	86 AEC 100 MW: Specific investment
711	electrolysis	lifetime	25.0000	years	Danish Energy Agency, data_sheets_for_renewabl...	86 AEC 100 MW: Technical lifetime
712	fuel cell	FOM	5.0000	%/year	Danish Energy Agency, technology_data_for_el_a...	12 LT-PEMFC CHP: Fixed O&M
713	fuel cell	c_b	1.2500	50oC/100oC	Danish Energy Agency, technology_data_for_el_a...	12 LT-PEMFC CHP: Cb coefficient
714	fuel cell	efficiency	0.5000	per unit	Danish Energy Agency, technology_data_for_el_a...	12 LT-PEMFC CHP: Electricity efficiency, annu...
715	fuel cell	investment	950.0000	EUR/kW_e	Danish Energy Agency, technology_data_for_el_a...	12 LT-PEMFC CHP: Nominal investment
716	fuel cell	lifetime	10.0000	years	Danish Energy Agency, technology_data_for_el_a...	12 LT-PEMFC CHP: Technical lifetime
717	CCGT	CO2 intensity	0.1980	tCO2/MWh_th	Stoichiometric calculation with 50 GJ/t CH4	NaN
718	CCGT	fuel	23.8481	EUR/MWh_th	DIW (2013): Current and propsective costs of e...	Based on IEA 2011 data.
745	hydro	FOM	1.0000	%/year	DIW DataDoc http://hdl.handle.net/10419/80348	from old pypsa cost assumptions
746	hydro	efficiency	0.9000	per unit	DIW DataDoc http://hdl.handle.net/10419/80348	from old pypsa cost assumptions
747	hydro	investment	2208.1616	EUR/kWel	DIW DataDoc http://hdl.handle.net/10419/80348	from old pypsa cost assumptions
748	hydro	lifetime	80.0000	years	IEA2010	from old pypsa cost assumptions
760	hydrogen storage underground	FOM	0.0000	%/year	Danish Energy Agency, technology_data_catalogu...	151c Hydrogen Storage - Caverns: Fixed O&M
761	hydrogen storage underground	VOM	0.0000	EUR/MWh	Danish Energy Agency, technology_data_catalogu...	151c Hydrogen Storage - Caverns: Variable O&M
762	hydrogen storage underground	investment	1.5000	EUR/kWh	Danish Energy Agency, technology_data_catalogu...	151c Hydrogen Storage - Caverns: Specific inv...
763	hydrogen storage underground	lifetime	100.0000	years	Danish Energy Agency, technology_data_catalogu...	151c Hydrogen Storage - Caverns: Technical li...
814	nuclear	FOM	1.2700	%/year	Lazard's levelized cost of energy analysis - v...	U.S. specific costs including newly commission...
815	nuclear	VOM	3.6188	EUR/MWh_e	Lazard's levelized cost of energy analysis - v...	U.S. specific costs including newly commission...
816	nuclear	efficiency	0.3260	p.u.	Lazard's levelized cost of energy analysis - v...	Based on heat rate of 10.45 MMBtu/MWh_e and 3....
817	nuclear	fuel	3.3122	EUR/MWh_th	DIW (2013): Current and propsective costs of e...	Based on IEA 2011 data.
818	nuclear	investment	8769.6136	EUR/kW_e	Lazard's levelized cost of energy analysis - v...	U.S. specific costs including newly commission...
819	nuclear	lifetime	40.0000	years	Lazard's levelized cost of energy analysis - v...	NaN
858	solar	FOM	2.0400	%/year	Calculated. See 'further description'.	Mixed investment costs based on average of 50%...
859	solar	VOM	0.0100	EUR/MWhel	RES costs made up to fix curtailment order	from old pypsa cost assumptions
860	solar	investment	407.8706	EUR/kW_e	Calculated. See 'further description'.	Mixed investment costs based on average of 50%...
861	solar	lifetime	40.0000	years	Calculated. See 'further description'.	Mixed investment costs based on average of 50%...

	technology	value	unit	source	further description
parameter
FOM	electrolysis	4.0000	%/year	Danish Energy Agency, data_sheets_for_renewabl...	86 AEC 100 MW: Fixed O&M
efficiency	electrolysis	0.6532	per unit	Danish Energy Agency, data_sheets_for_renewabl...	86 AEC 100 MW: Hydrogen Output
efficiency-heat	electrolysis	0.1849	per unit	Danish Energy Agency, data_sheets_for_renewabl...	86 AEC 100 MW: - hereof recoverable for dist...
investment	electrolysis	384.9356	EUR/kW_e	Danish Energy Agency, data_sheets_for_renewabl...	86 AEC 100 MW: Specific investment
lifetime	electrolysis	25.0000	years	Danish Energy Agency, data_sheets_for_renewabl...	86 AEC 100 MW: Technical lifetime

Load	Demand_load_Banskobystrický	Demand_load_Bratislavský	Demand_load_Košický	Demand_load_Nitriansky	Demand_load_Prešovský	Demand_load_Trenčiansky	Demand_load_Trnavský	Demand_load_Žilinský
snapshot
2017-01-01 00:00:00	364.359299	429.586052	459.745013	395.570447	476.604188	336.578964	333.569851	405.838708
2017-01-01 01:00:00	363.843475	428.977887	459.094152	395.010438	475.929459	336.102469	333.097616	405.264162
2017-01-01 02:00:00	362.529793	427.429033	457.436561	393.584226	474.211084	334.888948	331.894944	403.800928
2017-01-01 03:00:00	365.231299	430.614155	460.845294	396.517143	477.744817	337.384479	334.368164	406.809978
2017-01-01 04:00:00	381.691731	450.021296	481.614908	414.387582	499.276067	352.589897	349.437641	425.144300
...	...	...	...	...	...	...	...	...
2017-12-31 19:00:00	460.693214	543.165440	581.298209	500.156360	602.614826	425.567964	421.763263	513.139474
2017-12-31 20:00:00	435.896276	513.929412	550.009674	473.235307	570.178918	402.661652	399.061740	485.519601
2017-12-31 21:00:00	405.975189	478.651922	512.255539	440.751170	531.040312	375.021879	371.669074	452.192237
2017-12-31 22:00:00	377.655121	445.262061	476.521552	410.005194	493.995935	348.861055	345.742136	420.648154
2017-12-31 23:00:00	372.171151	438.796364	469.601932	404.051465	486.822567	343.795207	340.721578	414.539877

Assignment 4: System Planning¶

First, a land eligibility analysis, based on the criteria listede belo, is performed as follows:¶

Secondly, another land eligibility has been performed as follows, but this time using historical weather data of ERA5 dataset & atlite library.¶

Building the Model:¶

Below, a PyPSA model is being created that minimizes total annual system costs with the specified characteristics, according to the assignment explanations:¶

1) Solving The Model without a limit on CO2 emissions:¶

2) Solving the model with a CO2 emission reduction of 100% (no emissions):¶

	GID_0	GID_1	COUNTRY	geometry
NAME_1
Banskobystrický	SVK	SVK.1_1	Slovakia	MULTIPOLYGON (((19.45671 48.09119, 19.45659 48...
Bratislavský	SVK	SVK.2_1	Slovakia	MULTIPOLYGON (((17.31282 48.45557, 17.31436 48...
Košický	SVK	SVK.3_1	Slovakia	MULTIPOLYGON (((22.10827 48.38136, 22.10702 48...
Nitriansky	SVK	SVK.4_1	Slovakia	MULTIPOLYGON (((18.83296 47.81751, 18.81611 47...
Prešovský	SVK	SVK.5_1	Slovakia	MULTIPOLYGON (((21.54594 48.79453, 21.54533 48...
Trenčiansky	SVK	SVK.6_1	Slovakia	MULTIPOLYGON (((17.41709 48.83292, 17.42872 48...
Trnavský	SVK	SVK.7_1	Slovakia	MULTIPOLYGON (((17.45999 48.14781, 17.46112 48...
Žilinský	SVK	SVK.8_1	Slovakia	MULTIPOLYGON (((19.02072 48.87508, 19.01839 48...

NAME_1	Banskobystrický	Bratislavský	Košický	Nitriansky	Prešovský	Trenčiansky	Trnavský	Žilinský
time
2017-12-31 09:00:00	0.028286	0.013923	0.074632	0.020029	0.142586	0.017215	0.014964	0.035910
2017-12-31 10:00:00	0.049915	0.022226	0.118199	0.042010	0.244426	0.039485	0.027253	0.065328
2017-12-31 11:00:00	0.156638	0.026005	0.158670	0.046129	0.304047	0.045949	0.020684	0.143634
2017-12-31 12:00:00	0.085264	0.050059	0.127709	0.042245	0.222590	0.027803	0.034511	0.067804
2017-12-31 13:00:00	0.050969	0.161593	0.068615	0.061506	0.163986	0.024984	0.110701	0.047229
2017-12-31 14:00:00	0.038435	0.237795	0.022509	0.117132	0.056720	0.030638	0.184676	0.016817
2017-12-31 15:00:00	0.012403	0.129751	0.005537	0.069491	0.010469	0.028245	0.108905	0.003444
2017-12-31 16:00:00	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
2017-12-31 17:00:00	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
2017-12-31 18:00:00	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
2017-12-31 19:00:00	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
2017-12-31 20:00:00	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
2017-12-31 21:00:00	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
2017-12-31 22:00:00	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
2017-12-31 23:00:00	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000

attribute	v_nom	x	y	carrier	v_mag_pu_set	v_mag_pu_min	v_mag_pu_max	control
Bus
Banskobystrický	1.0	19.424031	48.499321	AC	1.0	0.0	inf	PQ
Bratislavský	1.0	17.175805	48.331983	AC	1.0	0.0	inf	PQ
Košický	1.0	21.292701	48.676239	AC	1.0	0.0	inf	PQ
Nitriansky	1.0	18.336016	48.223509	AC	1.0	0.0	inf	PQ
Prešovský	1.0	21.219622	49.114946	AC	1.0	0.0	inf	PQ
Trenčiansky	1.0	18.215575	48.902964	AC	1.0	0.0	inf	PQ
Trnavský	1.0	17.645627	48.323494	AC	1.0	0.0	inf	PQ
Žilinský	1.0	19.196299	49.177544	AC	1.0	0.0	inf	PQ

attribute	bus0	bus1	efficiency	build_year	lifetime	p_nom	p_nom_extendable	p_nom_min	...	shut_down_cost	min_up_time	min_down_time	up_time_before	down_time_before	ramp_limit_up	ramp_limit_down	ramp_limit_start_up	ramp_limit_shut_down	p_nom_opt
Link
TL1_Bratislavský_Trnavský	Bratislavský	Trnavský	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL2_Trnavský_Bratislavský	Trnavský	Bratislavský	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL1_Trnavský_Nitriansky	Trnavský	Nitriansky	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL2_Nitriansky_Trnavský	Nitriansky	Trnavský	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL1_Trnavský_Trenčiansky	Trnavský	Trenčiansky	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL2_Trenčiansky_Trnavský	Trenčiansky	Trnavský	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL1_Nitriansky_Trenčiansky	Nitriansky	Trenčiansky	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL2_Trenčiansky_Nitriansky	Trenčiansky	Nitriansky	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL1_Nitriansky_Banskobystrický	Nitriansky	Banskobystrický	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL2_Banskobystrický_Nitriansky	Banskobystrický	Nitriansky	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL1_Banskobystrický_Trenčiansky	Banskobystrický	Trenčiansky	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL2_Trenčiansky_Banskobystrický	Trenčiansky	Banskobystrický	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL1_Žilinský_Trenčiansky	Žilinský	Trenčiansky	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL2_Trenčiansky_Žilinský	Trenčiansky	Žilinský	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL1_Banskobystrický_Žilinský	Banskobystrický	Žilinský	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL2_Žilinský_Banskobystrický	Žilinský	Banskobystrický	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL1_Prešovský_Žilinský	Prešovský	Žilinský	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL2_Žilinský_Prešovský	Žilinský	Prešovský	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL1_Prešovský_Banskobystrický	Prešovský	Banskobystrický	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL2_Banskobystrický_Prešovský	Banskobystrický	Prešovský	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL1_Košický_Banskobystrický	Košický	Banskobystrický	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL2_Banskobystrický_Košický	Banskobystrický	Košický	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL1_Košický_Prešovský	Košický	Prešovský	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0
TL2_Prešovský_Košický	Prešovský	Košický	1.0	0	inf	0.0	False	0.0	...	0.0	0	0	1	0	NaN	NaN	1.0	1.0	0.0

attribute	bus	control	p_nom	p_nom_extendable	p_nom_min	p_nom_max	p_min_pu	p_max_pu	p_set	...	state_of_charge_initial_per_period	state_of_charge_set	cyclic_state_of_charge	cyclic_state_of_charge_per_period	max_hours	efficiency_store	efficiency_dispatch	standing_loss	inflow	p_nom_opt
StorageUnit
Battery_storage_Banskobystrický	Banskobystrický	PQ	0.0	True	0.0	inf	-1.0	1.0	0.0	...	False	NaN	True	True	6.0	0.96	0.96	0.0	0.0	0.0
Battery_storage_Bratislavský	Bratislavský	PQ	0.0	True	0.0	inf	-1.0	1.0	0.0	...	False	NaN	True	True	6.0	0.96	0.96	0.0	0.0	0.0
Battery_storage_Košický	Košický	PQ	0.0	True	0.0	inf	-1.0	1.0	0.0	...	False	NaN	True	True	6.0	0.96	0.96	0.0	0.0	0.0
Battery_storage_Nitriansky	Nitriansky	PQ	0.0	True	0.0	inf	-1.0	1.0	0.0	...	False	NaN	True	True	6.0	0.96	0.96	0.0	0.0	0.0
Battery_storage_Prešovský	Prešovský	PQ	0.0	True	0.0	inf	-1.0	1.0	0.0	...	False	NaN	True	True	6.0	0.96	0.96	0.0	0.0	0.0
Battery_storage_Trenčiansky	Trenčiansky	PQ	0.0	True	0.0	inf	-1.0	1.0	0.0	...	False	NaN	True	True	6.0	0.96	0.96	0.0	0.0	0.0
Battery_storage_Trnavský	Trnavský	PQ	0.0	True	0.0	inf	-1.0	1.0	0.0	...	False	NaN	True	True	6.0	0.96	0.96	0.0	0.0	0.0
Battery_storage_Žilinský	Žilinský	PQ	0.0	True	0.0	inf	-1.0	1.0	0.0	...	False	NaN	True	True	6.0	0.96	0.96	0.0	0.0	0.0

Generator	Coal_Košický	Coal_Trenčiansky	Gas_Bratislavský	Gas_Nitriansky	Gas_Trnavský	Hydro_Banskobystrický	Hydro_Bratislavský	Hydro_Košický	Hydro_Trenčiansky	Hydro_Trnavský	...	Solar_Žilinský	Wind_Žilinský	Banskobystrický load shedding	Bratislavský load shedding	Košický load shedding	Nitriansky load shedding	Prešovský load shedding	Trenčiansky load shedding	Trnavský load shedding	Žilinský load shedding
snapshot
2017-01-01 00:00:00	-0.0	-0.0	-0.0	-0.0	-0.0	14.381279	14.923516	18.437215	6.369863	247.155251	...	-0.0	271.806237	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0
2017-01-01 01:00:00	-0.0	-0.0	-0.0	-0.0	-0.0	14.381279	14.923516	18.437215	6.369863	247.155251	...	-0.0	295.034377	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0
2017-01-01 02:00:00	-0.0	-0.0	-0.0	-0.0	-0.0	14.381279	14.923516	18.437215	6.369863	247.155251	...	-0.0	313.706895	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0
2017-01-01 03:00:00	-0.0	-0.0	-0.0	-0.0	-0.0	14.381279	14.923516	18.437215	6.369863	247.155251	...	-0.0	362.456009	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0
2017-01-01 04:00:00	-0.0	-0.0	-0.0	-0.0	-0.0	14.381279	14.923516	18.437215	6.369863	247.155251	...	-0.0	423.834183	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
2017-12-31 19:00:00	-0.0	-0.0	-0.0	-0.0	-0.0	14.381279	-0.000000	18.437215	-0.000000	-0.000000	...	-0.0	-0.000000	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0
2017-12-31 20:00:00	-0.0	-0.0	-0.0	-0.0	-0.0	14.381279	14.923516	18.437215	6.369863	-0.000000	...	-0.0	485.519601	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0
2017-12-31 21:00:00	-0.0	-0.0	-0.0	-0.0	-0.0	14.381279	14.923516	18.437215	6.369863	-0.000000	...	-0.0	-0.000000	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0
2017-12-31 22:00:00	-0.0	-0.0	-0.0	-0.0	-0.0	14.381279	14.923516	18.437215	6.369863	-0.000000	...	-0.0	-0.000000	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0
2017-12-31 23:00:00	-0.0	-0.0	-0.0	-0.0	-0.0	14.381279	-0.000000	18.437215	6.369863	-0.000000	...	-0.0	414.539877	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0	-0.0